<?php

/**
 * @package JAMA
 *
 *    Class to obtain eigenvalues and eigenvectors of a real matrix.
 *
 *    If A is symmetric, then A = V*D*V' where the eigenvalue matrix D
 *    is diagonal and the eigenvector matrix V is orthogonal (i.e.
 *    A = V.times(D.times(V.transpose())) and V.times(V.transpose())
 *    equals the identity matrix).
 *
 *    If A is not symmetric, then the eigenvalue matrix D is block diagonal
 *    with the real eigenvalues in 1-by-1 blocks and any complex eigenvalues,
 *    lambda + i*mu, in 2-by-2 blocks, [lambda, mu; -mu, lambda].  The
 *    columns of V represent the eigenvectors in the sense that A*V = V*D,
 *    i.e. A.times(V) equals V.times(D).  The matrix V may be badly
 *    conditioned, or even singular, so the validity of the equation
 *    A = V*D*inverse(V) depends upon V.cond().
 *
 * @author  Paul Meagher
 * @license PHP v3.0
 * @version 1.1
 */
class EigenvalueDecomposition {

	/**
	 *    Row and column dimension (square matrix).
	 *
	 * @var int
	 */
	private $n;

	/**
	 *    Internal symmetry flag.
	 *
	 * @var int
	 */
	private $issymmetric;

	/**
	 *    Arrays for internal storage of eigenvalues.
	 *
	 * @var array
	 */
	private $d = array();
	private $e = array();

	/**
	 *    Array for internal storage of eigenvectors.
	 *
	 * @var array
	 */
	private $V = array();

	/**
	 *    Array for internal storage of nonsymmetric Hessenberg form.
	 *
	 * @var array
	 */
	private $H = array();

	/**
	 *    Working storage for nonsymmetric algorithm.
	 *
	 * @var array
	 */
	private $ort;

	/**
	 *    Used for complex scalar division.
	 *
	 * @var float
	 */
	private $cdivr;
	private $cdivi;


	/**
	 *    Symmetric Householder reduction to tridiagonal form.
	 *
	 * @access private
	 */
	private function tred2() {
		//  This is derived from the Algol procedures tred2 by
		//  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
		//  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
		//  Fortran subroutine in EISPACK.
		$this->d = $this->V[$this->n - 1];
		// Householder reduction to tridiagonal form.
		for ($i = $this->n - 1; $i > 0; --$i)
		{
			$i_ = $i - 1;
			// Scale to avoid under/overflow.
			$h = $scale = 0.0;
			$scale += array_sum(array_map(abs, $this->d));
			if ($scale == 0.0)
			{
				$this->e[$i] = $this->d[$i_];
				$this->d = array_slice($this->V[$i_], 0, $i_);
				for ($j = 0; $j < $i; ++$j)
				{
					$this->V[$j][$i] = $this->V[$i][$j] = 0.0;
				}
			}
			else
			{
				// Generate Householder vector.
				for ($k = 0; $k < $i; ++$k)
				{
					$this->d[$k] /= $scale;
					$h += pow($this->d[$k], 2);
				}
				$f = $this->d[$i_];
				$g = sqrt($h);
				if ($f > 0)
				{
					$g = -$g;
				}
				$this->e[$i] = $scale * $g;
				$h = $h - $f * $g;
				$this->d[$i_] = $f - $g;
				for ($j = 0; $j < $i; ++$j)
				{
					$this->e[$j] = 0.0;
				}
				// Apply similarity transformation to remaining columns.
				for ($j = 0; $j < $i; ++$j)
				{
					$f = $this->d[$j];
					$this->V[$j][$i] = $f;
					$g = $this->e[$j] + $this->V[$j][$j] * $f;
					for ($k = $j + 1; $k <= $i_; ++$k)
					{
						$g += $this->V[$k][$j] * $this->d[$k];
						$this->e[$k] += $this->V[$k][$j] * $f;
					}
					$this->e[$j] = $g;
				}
				$f = 0.0;
				for ($j = 0; $j < $i; ++$j)
				{
					$this->e[$j] /= $h;
					$f += $this->e[$j] * $this->d[$j];
				}
				$hh = $f / (2 * $h);
				for ($j = 0; $j < $i; ++$j)
				{
					$this->e[$j] -= $hh * $this->d[$j];
				}
				for ($j = 0; $j < $i; ++$j)
				{
					$f = $this->d[$j];
					$g = $this->e[$j];
					for ($k = $j; $k <= $i_; ++$k)
					{
						$this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
					}
					$this->d[$j] = $this->V[$i - 1][$j];
					$this->V[$i][$j] = 0.0;
				}
			}
			$this->d[$i] = $h;
		}

		// Accumulate transformations.
		for ($i = 0; $i < $this->n - 1; ++$i)
		{
			$this->V[$this->n - 1][$i] = $this->V[$i][$i];
			$this->V[$i][$i] = 1.0;
			$h = $this->d[$i + 1];
			if ($h != 0.0)
			{
				for ($k = 0; $k <= $i; ++$k)
				{
					$this->d[$k] = $this->V[$k][$i + 1] / $h;
				}
				for ($j = 0; $j <= $i; ++$j)
				{
					$g = 0.0;
					for ($k = 0; $k <= $i; ++$k)
					{
						$g += $this->V[$k][$i + 1] * $this->V[$k][$j];
					}
					for ($k = 0; $k <= $i; ++$k)
					{
						$this->V[$k][$j] -= $g * $this->d[$k];
					}
				}
			}
			for ($k = 0; $k <= $i; ++$k)
			{
				$this->V[$k][$i + 1] = 0.0;
			}
		}

		$this->d = $this->V[$this->n - 1];
		$this->V[$this->n - 1] = array_fill(0, $j, 0.0);
		$this->V[$this->n - 1][$this->n - 1] = 1.0;
		$this->e[0] = 0.0;
	}


	/**
	 *    Symmetric tridiagonal QL algorithm.
	 *
	 *    This is derived from the Algol procedures tql2, by
	 *    Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
	 *    Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
	 *    Fortran subroutine in EISPACK.
	 *
	 * @access private
	 */
	private function tql2() {
		for ($i = 1; $i < $this->n; ++$i)
		{
			$this->e[$i - 1] = $this->e[$i];
		}
		$this->e[$this->n - 1] = 0.0;
		$f = 0.0;
		$tst1 = 0.0;
		$eps = pow(2.0, -52.0);

		for ($l = 0; $l < $this->n; ++$l)
		{
			// Find small subdiagonal element
			$tst1 = max($tst1, abs($this->d[$l]) + abs($this->e[$l]));
			$m = $l;
			while ($m < $this->n)
			{
				if (abs($this->e[$m]) <= $eps * $tst1)
					break;
				++$m;
			}
			// If m == l, $this->d[l] is an eigenvalue,
			// otherwise, iterate.
			if ($m > $l)
			{
				$iter = 0;
				do
				{
					// Could check iteration count here.
					$iter += 1;
					// Compute implicit shift
					$g = $this->d[$l];
					$p = ($this->d[$l + 1] - $g) / (2.0 * $this->e[$l]);
					$r = hypo($p, 1.0);
					if ($p < 0)
						$r *= -1;
					$this->d[$l] = $this->e[$l] / ($p + $r);
					$this->d[$l + 1] = $this->e[$l] * ($p + $r);
					$dl1 = $this->d[$l + 1];
					$h = $g - $this->d[$l];
					for ($i = $l + 2; $i < $this->n; ++$i)
						$this->d[$i] -= $h;
					$f += $h;
					// Implicit QL transformation.
					$p = $this->d[$m];
					$c = 1.0;
					$c2 = $c3 = $c;
					$el1 = $this->e[$l + 1];
					$s = $s2 = 0.0;
					for ($i = $m - 1; $i >= $l; --$i)
					{
						$c3 = $c2;
						$c2 = $c;
						$s2 = $s;
						$g = $c * $this->e[$i];
						$h = $c * $p;
						$r = hypo($p, $this->e[$i]);
						$this->e[$i + 1] = $s * $r;
						$s = $this->e[$i] / $r;
						$c = $p / $r;
						$p = $c * $this->d[$i] - $s * $g;
						$this->d[$i + 1] = $h + $s * ($c * $g + $s * $this->d[$i]);
						// Accumulate transformation.
						for ($k = 0; $k < $this->n; ++$k)
						{
							$h = $this->V[$k][$i + 1];
							$this->V[$k][$i + 1] = $s * $this->V[$k][$i] + $c * $h;
							$this->V[$k][$i] = $c * $this->V[$k][$i] - $s * $h;
						}
					}
					$p = -$s * $s2 * $c3 * $el1 * $this->e[$l] / $dl1;
					$this->e[$l] = $s * $p;
					$this->d[$l] = $c * $p;
					// Check for convergence.
				} while (abs($this->e[$l]) > $eps * $tst1);
			}
			$this->d[$l] = $this->d[$l] + $f;
			$this->e[$l] = 0.0;
		}

		// Sort eigenvalues and corresponding vectors.
		for ($i = 0; $i < $this->n - 1; ++$i)
		{
			$k = $i;
			$p = $this->d[$i];
			for ($j = $i + 1; $j < $this->n; ++$j)
			{
				if ($this->d[$j] < $p)
				{
					$k = $j;
					$p = $this->d[$j];
				}
			}
			if ($k != $i)
			{
				$this->d[$k] = $this->d[$i];
				$this->d[$i] = $p;
				for ($j = 0; $j < $this->n; ++$j)
				{
					$p = $this->V[$j][$i];
					$this->V[$j][$i] = $this->V[$j][$k];
					$this->V[$j][$k] = $p;
				}
			}
		}
	}


	/**
	 *    Nonsymmetric reduction to Hessenberg form.
	 *
	 *    This is derived from the Algol procedures orthes and ortran,
	 *    by Martin and Wilkinson, Handbook for Auto. Comp.,
	 *    Vol.ii-Linear Algebra, and the corresponding
	 *    Fortran subroutines in EISPACK.
	 *
	 * @access private
	 */
	private function orthes() {
		$low = 0;
		$high = $this->n - 1;

		for ($m = $low + 1; $m <= $high - 1; ++$m)
		{
			// Scale column.
			$scale = 0.0;
			for ($i = $m; $i <= $high; ++$i)
			{
				$scale = $scale + abs($this->H[$i][$m - 1]);
			}
			if ($scale != 0.0)
			{
				// Compute Householder transformation.
				$h = 0.0;
				for ($i = $high; $i >= $m; --$i)
				{
					$this->ort[$i] = $this->H[$i][$m - 1] / $scale;
					$h += $this->ort[$i] * $this->ort[$i];
				}
				$g = sqrt($h);
				if ($this->ort[$m] > 0)
				{
					$g *= -1;
				}
				$h -= $this->ort[$m] * $g;
				$this->ort[$m] -= $g;
				// Apply Householder similarity transformation
				// H = (I -u * u' / h) * H * (I -u * u') / h)
				for ($j = $m; $j < $this->n; ++$j)
				{
					$f = 0.0;
					for ($i = $high; $i >= $m; --$i)
					{
						$f += $this->ort[$i] * $this->H[$i][$j];
					}
					$f /= $h;
					for ($i = $m; $i <= $high; ++$i)
					{
						$this->H[$i][$j] -= $f * $this->ort[$i];
					}
				}
				for ($i = 0; $i <= $high; ++$i)
				{
					$f = 0.0;
					for ($j = $high; $j >= $m; --$j)
					{
						$f += $this->ort[$j] * $this->H[$i][$j];
					}
					$f = $f / $h;
					for ($j = $m; $j <= $high; ++$j)
					{
						$this->H[$i][$j] -= $f * $this->ort[$j];
					}
				}
				$this->ort[$m] = $scale * $this->ort[$m];
				$this->H[$m][$m - 1] = $scale * $g;
			}
		}

		// Accumulate transformations (Algol's ortran).
		for ($i = 0; $i < $this->n; ++$i)
		{
			for ($j = 0; $j < $this->n; ++$j)
			{
				$this->V[$i][$j] = ($i == $j ? 1.0 : 0.0);
			}
		}
		for ($m = $high - 1; $m >= $low + 1; --$m)
		{
			if ($this->H[$m][$m - 1] != 0.0)
			{
				for ($i = $m + 1; $i <= $high; ++$i)
				{
					$this->ort[$i] = $this->H[$i][$m - 1];
				}
				for ($j = $m; $j <= $high; ++$j)
				{
					$g = 0.0;
					for ($i = $m; $i <= $high; ++$i)
					{
						$g += $this->ort[$i] * $this->V[$i][$j];
					}
					// Double division avoids possible underflow
					$g = ($g / $this->ort[$m]) / $this->H[$m][$m - 1];
					for ($i = $m; $i <= $high; ++$i)
					{
						$this->V[$i][$j] += $g * $this->ort[$i];
					}
				}
			}
		}
	}


	/**
	 *    Performs complex division.
	 *
	 * @access private
	 */
	private function cdiv($xr, $xi, $yr, $yi) {
		if (abs($yr) > abs($yi))
		{
			$r = $yi / $yr;
			$d = $yr + $r * $yi;
			$this->cdivr = ($xr + $r * $xi) / $d;
			$this->cdivi = ($xi - $r * $xr) / $d;
		}
		else
		{
			$r = $yr / $yi;
			$d = $yi + $r * $yr;
			$this->cdivr = ($r * $xr + $xi) / $d;
			$this->cdivi = ($r * $xi - $xr) / $d;
		}
	}


	/**
	 *    Nonsymmetric reduction from Hessenberg to real Schur form.
	 *
	 *    Code is derived from the Algol procedure hqr2,
	 *    by Martin and Wilkinson, Handbook for Auto. Comp.,
	 *    Vol.ii-Linear Algebra, and the corresponding
	 *    Fortran subroutine in EISPACK.
	 *
	 * @access private
	 */
	private function hqr2() {
		//  Initialize
		$nn = $this->n;
		$n = $nn - 1;
		$low = 0;
		$high = $nn - 1;
		$eps = pow(2.0, -52.0);
		$exshift = 0.0;
		$p = $q = $r = $s = $z = 0;
		// Store roots isolated by balanc and compute matrix norm
		$norm = 0.0;

		for ($i = 0; $i < $nn; ++$i)
		{
			if (($i < $low) OR ($i > $high))
			{
				$this->d[$i] = $this->H[$i][$i];
				$this->e[$i] = 0.0;
			}
			for ($j = max($i - 1, 0); $j < $nn; ++$j)
			{
				$norm = $norm + abs($this->H[$i][$j]);
			}
		}

		// Outer loop over eigenvalue index
		$iter = 0;
		while ($n >= $low)
		{
			// Look for single small sub-diagonal element
			$l = $n;
			while ($l > $low)
			{
				$s = abs($this->H[$l - 1][$l - 1]) + abs($this->H[$l][$l]);
				if ($s == 0.0)
				{
					$s = $norm;
				}
				if (abs($this->H[$l][$l - 1]) < $eps * $s)
				{
					break;
				}
				--$l;
			}
			// Check for convergence
			// One root found
			if ($l == $n)
			{
				$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
				$this->d[$n] = $this->H[$n][$n];
				$this->e[$n] = 0.0;
				--$n;
				$iter = 0;
				// Two roots found
			}
			else if ($l == $n - 1)
			{
				$w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
				$p = ($this->H[$n - 1][$n - 1] - $this->H[$n][$n]) / 2.0;
				$q = $p * $p + $w;
				$z = sqrt(abs($q));
				$this->H[$n][$n] = $this->H[$n][$n] + $exshift;
				$this->H[$n - 1][$n - 1] = $this->H[$n - 1][$n - 1] + $exshift;
				$x = $this->H[$n][$n];
				// Real pair
				if ($q >= 0)
				{
					if ($p >= 0)
					{
						$z = $p + $z;
					}
					else
					{
						$z = $p - $z;
					}
					$this->d[$n - 1] = $x + $z;
					$this->d[$n] = $this->d[$n - 1];
					if ($z != 0.0)
					{
						$this->d[$n] = $x - $w / $z;
					}
					$this->e[$n - 1] = 0.0;
					$this->e[$n] = 0.0;
					$x = $this->H[$n][$n - 1];
					$s = abs($x) + abs($z);
					$p = $x / $s;
					$q = $z / $s;
					$r = sqrt($p * $p + $q * $q);
					$p = $p / $r;
					$q = $q / $r;
					// Row modification
					for ($j = $n - 1; $j < $nn; ++$j)
					{
						$z = $this->H[$n - 1][$j];
						$this->H[$n - 1][$j] = $q * $z + $p * $this->H[$n][$j];
						$this->H[$n][$j] = $q * $this->H[$n][$j] - $p * $z;
					}
					// Column modification
					for ($i = 0; $i <= n; ++$i)
					{
						$z = $this->H[$i][$n - 1];
						$this->H[$i][$n - 1] = $q * $z + $p * $this->H[$i][$n];
						$this->H[$i][$n] = $q * $this->H[$i][$n] - $p * $z;
					}
					// Accumulate transformations
					for ($i = $low; $i <= $high; ++$i)
					{
						$z = $this->V[$i][$n - 1];
						$this->V[$i][$n - 1] = $q * $z + $p * $this->V[$i][$n];
						$this->V[$i][$n] = $q * $this->V[$i][$n] - $p * $z;
					}
					// Complex pair
				}
				else
				{
					$this->d[$n - 1] = $x + $p;
					$this->d[$n] = $x + $p;
					$this->e[$n - 1] = $z;
					$this->e[$n] = -$z;
				}
				$n = $n - 2;
				$iter = 0;
				// No convergence yet
			}
			else
			{
				// Form shift
				$x = $this->H[$n][$n];
				$y = 0.0;
				$w = 0.0;
				if ($l < $n)
				{
					$y = $this->H[$n - 1][$n - 1];
					$w = $this->H[$n][$n - 1] * $this->H[$n - 1][$n];
				}
				// Wilkinson's original ad hoc shift
				if ($iter == 10)
				{
					$exshift += $x;
					for ($i = $low; $i <= $n; ++$i)
					{
						$this->H[$i][$i] -= $x;
					}
					$s = abs($this->H[$n][$n - 1]) + abs($this->H[$n - 1][$n - 2]);
					$x = $y = 0.75 * $s;
					$w = -0.4375 * $s * $s;
				}
				// MATLAB's new ad hoc shift
				if ($iter == 30)
				{
					$s = ($y - $x) / 2.0;
					$s = $s * $s + $w;
					if ($s > 0)
					{
						$s = sqrt($s);
						if ($y < $x)
						{
							$s = -$s;
						}
						$s = $x - $w / (($y - $x) / 2.0 + $s);
						for ($i = $low; $i <= $n; ++$i)
						{
							$this->H[$i][$i] -= $s;
						}
						$exshift += $s;
						$x = $y = $w = 0.964;
					}
				}
				// Could check iteration count here.
				$iter = $iter + 1;
				// Look for two consecutive small sub-diagonal elements
				$m = $n - 2;
				while ($m >= $l)
				{
					$z = $this->H[$m][$m];
					$r = $x - $z;
					$s = $y - $z;
					$p = ($r * $s - $w) / $this->H[$m + 1][$m] + $this->H[$m][$m + 1];
					$q = $this->H[$m + 1][$m + 1] - $z - $r - $s;
					$r = $this->H[$m + 2][$m + 1];
					$s = abs($p) + abs($q) + abs($r);
					$p = $p / $s;
					$q = $q / $s;
					$r = $r / $s;
					if ($m == $l)
					{
						break;
					}
					if (abs($this->H[$m][$m - 1]) * (abs($q) + abs($r)) < $eps * (abs($p) * (abs($this->H[$m - 1][$m - 1]) + abs($z) + abs($this->H[$m + 1][$m + 1]))))
					{
						break;
					}
					--$m;
				}
				for ($i = $m + 2; $i <= $n; ++$i)
				{
					$this->H[$i][$i - 2] = 0.0;
					if ($i > $m + 2)
					{
						$this->H[$i][$i - 3] = 0.0;
					}
				}
				// Double QR step involving rows l:n and columns m:n
				for ($k = $m; $k <= $n - 1; ++$k)
				{
					$notlast = ($k != $n - 1);
					if ($k != $m)
					{
						$p = $this->H[$k][$k - 1];
						$q = $this->H[$k + 1][$k - 1];
						$r = ($notlast ? $this->H[$k + 2][$k - 1] : 0.0);
						$x = abs($p) + abs($q) + abs($r);
						if ($x != 0.0)
						{
							$p = $p / $x;
							$q = $q / $x;
							$r = $r / $x;
						}
					}
					if ($x == 0.0)
					{
						break;
					}
					$s = sqrt($p * $p + $q * $q + $r * $r);
					if ($p < 0)
					{
						$s = -$s;
					}
					if ($s != 0)
					{
						if ($k != $m)
						{
							$this->H[$k][$k - 1] = -$s * $x;
						}
						elseif ($l != $m)
						{
							$this->H[$k][$k - 1] = -$this->H[$k][$k - 1];
						}
						$p = $p + $s;
						$x = $p / $s;
						$y = $q / $s;
						$z = $r / $s;
						$q = $q / $p;
						$r = $r / $p;
						// Row modification
						for ($j = $k; $j < $nn; ++$j)
						{
							$p = $this->H[$k][$j] + $q * $this->H[$k + 1][$j];
							if ($notlast)
							{
								$p = $p + $r * $this->H[$k + 2][$j];
								$this->H[$k + 2][$j] = $this->H[$k + 2][$j] - $p * $z;
							}
							$this->H[$k][$j] = $this->H[$k][$j] - $p * $x;
							$this->H[$k + 1][$j] = $this->H[$k + 1][$j] - $p * $y;
						}
						// Column modification
						for ($i = 0; $i <= min($n, $k + 3); ++$i)
						{
							$p = $x * $this->H[$i][$k] + $y * $this->H[$i][$k + 1];
							if ($notlast)
							{
								$p = $p + $z * $this->H[$i][$k + 2];
								$this->H[$i][$k + 2] = $this->H[$i][$k + 2] - $p * $r;
							}
							$this->H[$i][$k] = $this->H[$i][$k] - $p;
							$this->H[$i][$k + 1] = $this->H[$i][$k + 1] - $p * $q;
						}
						// Accumulate transformations
						for ($i = $low; $i <= $high; ++$i)
						{
							$p = $x * $this->V[$i][$k] + $y * $this->V[$i][$k + 1];
							if ($notlast)
							{
								$p = $p + $z * $this->V[$i][$k + 2];
								$this->V[$i][$k + 2] = $this->V[$i][$k + 2] - $p * $r;
							}
							$this->V[$i][$k] = $this->V[$i][$k] - $p;
							$this->V[$i][$k + 1] = $this->V[$i][$k + 1] - $p * $q;
						}
					} // ($s != 0)
				} // k loop
			} // check convergence
		} // while ($n >= $low)

		// Backsubstitute to find vectors of upper triangular form
		if ($norm == 0.0)
		{
			return;
		}

		for ($n = $nn - 1; $n >= 0; --$n)
		{
			$p = $this->d[$n];
			$q = $this->e[$n];
			// Real vector
			if ($q == 0)
			{
				$l = $n;
				$this->H[$n][$n] = 1.0;
				for ($i = $n - 1; $i >= 0; --$i)
				{
					$w = $this->H[$i][$i] - $p;
					$r = 0.0;
					for ($j = $l; $j <= $n; ++$j)
					{
						$r = $r + $this->H[$i][$j] * $this->H[$j][$n];
					}
					if ($this->e[$i] < 0.0)
					{
						$z = $w;
						$s = $r;
					}
					else
					{
						$l = $i;
						if ($this->e[$i] == 0.0)
						{
							if ($w != 0.0)
							{
								$this->H[$i][$n] = -$r / $w;
							}
							else
							{
								$this->H[$i][$n] = -$r / ($eps * $norm);
							}
							// Solve real equations
						}
						else
						{
							$x = $this->H[$i][$i + 1];
							$y = $this->H[$i + 1][$i];
							$q = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i];
							$t = ($x * $s - $z * $r) / $q;
							$this->H[$i][$n] = $t;
							if (abs($x) > abs($z))
							{
								$this->H[$i + 1][$n] = (-$r - $w * $t) / $x;
							}
							else
							{
								$this->H[$i + 1][$n] = (-$s - $y * $t) / $z;
							}
						}
						// Overflow control
						$t = abs($this->H[$i][$n]);
						if (($eps * $t) * $t > 1)
						{
							for ($j = $i; $j <= $n; ++$j)
							{
								$this->H[$j][$n] = $this->H[$j][$n] / $t;
							}
						}
					}
				}
				// Complex vector
			}
			else if ($q < 0)
			{
				$l = $n - 1;
				// Last vector component imaginary so matrix is triangular
				if (abs($this->H[$n][$n - 1]) > abs($this->H[$n - 1][$n]))
				{
					$this->H[$n - 1][$n - 1] = $q / $this->H[$n][$n - 1];
					$this->H[$n - 1][$n] = -($this->H[$n][$n] - $p) / $this->H[$n][$n - 1];
				}
				else
				{
					$this->cdiv(0.0, -$this->H[$n - 1][$n], $this->H[$n - 1][$n - 1] - $p, $q);
					$this->H[$n - 1][$n - 1] = $this->cdivr;
					$this->H[$n - 1][$n] = $this->cdivi;
				}
				$this->H[$n][$n - 1] = 0.0;
				$this->H[$n][$n] = 1.0;
				for ($i = $n - 2; $i >= 0; --$i)
				{
					// double ra,sa,vr,vi;
					$ra = 0.0;
					$sa = 0.0;
					for ($j = $l; $j <= $n; ++$j)
					{
						$ra = $ra + $this->H[$i][$j] * $this->H[$j][$n - 1];
						$sa = $sa + $this->H[$i][$j] * $this->H[$j][$n];
					}
					$w = $this->H[$i][$i] - $p;
					if ($this->e[$i] < 0.0)
					{
						$z = $w;
						$r = $ra;
						$s = $sa;
					}
					else
					{
						$l = $i;
						if ($this->e[$i] == 0)
						{
							$this->cdiv(-$ra, -$sa, $w, $q);
							$this->H[$i][$n - 1] = $this->cdivr;
							$this->H[$i][$n] = $this->cdivi;
						}
						else
						{
							// Solve complex equations
							$x = $this->H[$i][$i + 1];
							$y = $this->H[$i + 1][$i];
							$vr = ($this->d[$i] - $p) * ($this->d[$i] - $p) + $this->e[$i] * $this->e[$i] - $q * $q;
							$vi = ($this->d[$i] - $p) * 2.0 * $q;
							if ($vr == 0.0 & $vi == 0.0)
							{
								$vr = $eps * $norm * (abs($w) + abs($q) + abs($x) + abs($y) + abs($z));
							}
							$this->cdiv($x * $r - $z * $ra + $q * $sa, $x * $s - $z * $sa - $q * $ra, $vr, $vi);
							$this->H[$i][$n - 1] = $this->cdivr;
							$this->H[$i][$n] = $this->cdivi;
							if (abs($x) > (abs($z) + abs($q)))
							{
								$this->H[$i + 1][$n - 1] = (-$ra - $w * $this->H[$i][$n - 1] + $q * $this->H[$i][$n]) / $x;
								$this->H[$i + 1][$n] = (-$sa - $w * $this->H[$i][$n] - $q * $this->H[$i][$n - 1]) / $x;
							}
							else
							{
								$this->cdiv(-$r - $y * $this->H[$i][$n - 1], -$s - $y * $this->H[$i][$n], $z, $q);
								$this->H[$i + 1][$n - 1] = $this->cdivr;
								$this->H[$i + 1][$n] = $this->cdivi;
							}
						}
						// Overflow control
						$t = max(abs($this->H[$i][$n - 1]), abs($this->H[$i][$n]));
						if (($eps * $t) * $t > 1)
						{
							for ($j = $i; $j <= $n; ++$j)
							{
								$this->H[$j][$n - 1] = $this->H[$j][$n - 1] / $t;
								$this->H[$j][$n] = $this->H[$j][$n] / $t;
							}
						}
					} // end else
				} // end for
			} // end else for complex case
		} // end for

		// Vectors of isolated roots
		for ($i = 0; $i < $nn; ++$i)
		{
			if ($i < $low | $i > $high)
			{
				for ($j = $i; $j < $nn; ++$j)
				{
					$this->V[$i][$j] = $this->H[$i][$j];
				}
			}
		}

		// Back transformation to get eigenvectors of original matrix
		for ($j = $nn - 1; $j >= $low; --$j)
		{
			for ($i = $low; $i <= $high; ++$i)
			{
				$z = 0.0;
				for ($k = $low; $k <= min($j, $high); ++$k)
				{
					$z = $z + $this->V[$i][$k] * $this->H[$k][$j];
				}
				$this->V[$i][$j] = $z;
			}
		}
	} // end hqr2


	/**
	 *    Constructor: Check for symmetry, then construct the eigenvalue decomposition
	 *
	 * @access public
	 *
	 * @param A  Square matrix
	 *
	 * @return Structure to access D and V.
	 */
	public function __construct($Arg) {
		$this->A = $Arg->getArray();
		$this->n = $Arg->getColumnDimension();

		$issymmetric = true;
		for ($j = 0; ($j < $this->n) & $issymmetric; ++$j)
		{
			for ($i = 0; ($i < $this->n) & $issymmetric; ++$i)
			{
				$issymmetric = ($this->A[$i][$j] == $this->A[$j][$i]);
			}
		}

		if ($issymmetric)
		{
			$this->V = $this->A;
			// Tridiagonalize.
			$this->tred2();
			// Diagonalize.
			$this->tql2();
		}
		else
		{
			$this->H = $this->A;
			$this->ort = array();
			// Reduce to Hessenberg form.
			$this->orthes();
			// Reduce Hessenberg to real Schur form.
			$this->hqr2();
		}
	}


	/**
	 *    Return the eigenvector matrix
	 *
	 * @access public
	 * @return V
	 */
	public function getV() {
		return new Matrix($this->V, $this->n, $this->n);
	}


	/**
	 *    Return the real parts of the eigenvalues
	 *
	 * @access public
	 * @return real(diag(D))
	 */
	public function getRealEigenvalues() {
		return $this->d;
	}


	/**
	 *    Return the imaginary parts of the eigenvalues
	 *
	 * @access public
	 * @return imag(diag(D))
	 */
	public function getImagEigenvalues() {
		return $this->e;
	}


	/**
	 *    Return the block diagonal eigenvalue matrix
	 *
	 * @access public
	 * @return D
	 */
	public function getD() {
		for ($i = 0; $i < $this->n; ++$i)
		{
			$D[$i] = array_fill(0, $this->n, 0.0);
			$D[$i][$i] = $this->d[$i];
			if ($this->e[$i] == 0)
			{
				continue;
			}
			$o = ($this->e[$i] > 0) ? $i + 1 : $i - 1;
			$D[$i][$o] = $this->e[$i];
		}

		return new Matrix($D);
	}

} //	class EigenvalueDecomposition
